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ABSTRACT 

We model gas cooling in high-resolution N-body simulations in order to investigate the 
formation of the first generation of stars. We follow a region of a Lambda Cold Dark 
Matter (ACDM) universe especially selected to contain a rich cluster by the present 
day. The properties of the dark haloes that form in these sub-solar mass resolution 
simulations are presented in a companion paper by Gao et al. The first gas clouds 
able to cool by molecular hydrogen line emission collapse at extremely high redshift, 
z » 47, when the mass of the dark halo is 2.4x 10 5 h -1 Mo. By z « 30, a substantial 
population of haloes are capable of undergoing molecular hydrogen cooling although 
their ability to form stars is dependent on the efficiency of feedback processes such 
as dissociating Lyman- Werner radiation. The mass of the main halo grows extremely 
rapidly and, by z w 36, its virial temperature has reached 10 4 K, at which point gas 
cooling becomes dominated by more effective atomic processes. By z « 30, a small 
"group" of such potential galaxies will have formed unless prevented from doing so 
by feedback processes. By this redshift, massive ( £ 100M Q ) population III stars are 
able to ionise gas well beyond their own host halo and neighbouring HII regions can 
percolate to form an ionised superbubble. Such patches would be too widely separated 
to contribute significantly to reionisation at this time. The large number density of 
early cooling haloes in the pre-reionised universe raises the exciting prospect that this 
ultra-early generation of stars may be observable as gamma-ray bursts or supernovae. 

Key words: galaxies: haloes - galaxies: formation - methods: N-body simulations - 
cosmology: theory - cosmology: dark matter 



1 INTRODUCTION 

The first stars in the universe are believed to have formed 
from primordial metal free gas in haloes with virial temper- 
atures less than the ~ 10 4 K threshold for atomic hydrogen 
line cooling (see reviews by Bromm & Larson 2004; Ciardi & 
Ferrara 2005 and references therein). Primordial gas in these 
small haloes is instead able to cool by molecular hydrogen 
transitions, a less effective process (Saslaw & Zipoy 1967; 
Peebles & Dicke 1968). Molecular hydrogen is produced in 
a two step reaction chain in which free electrons left over 
from recombination serve as a catalyst via H~ production 
(when z ,$ 100). Metal free (population III) stars are the 
first potential producers of UV photons that can contribute 
to the reionisation process, and are the first producers of the 
metals required for the formation of population II stars. 



* Email: d.s.reed@durham.ac.uk 



Tegmark et al. (1997) developed analytic methods to 
model early baryonic collapse via FL cooling. The early 
stages of primordial star formation were later directly mod- 
elled in high-resolution hydrodynamic numerical simulations 
by Bromm, Coppi & Larson (1999, 2002) and Abel, Bryan & 
Norman (2002). Yoshida et al. (2003) further utilised simula- 
tions to develop a semi-analytic model based on the Tegmark 
et al. methods and included the effects of dynamical heating 
caused by the thermalisation of kinetic energy of infall into 
a deepening halo potential. This heating effect varies from 
halo to halo, depending on mass accretion rates, and has the 
potential to delay baryonic collapse. 

Determining the first epoch of baryonic cooling in high 
redshift haloes thus requires the use of numerical simulations 
to follow their growth. In this paper, we study the growth 
of structure in a region selected to be the site where a rich 
cluster forms by z = in a large cosmological simulation. We 
repeatedly "resimulate" this region at increasing resolution 



© 2003 RAS 



2 Reed et al. 



and redshift. In this way, we are able to model the site of 
an extremely early generation of star formation. The sim- 
ulations employed here are those described in a companion 
paper (Gao et al. 2005; G05 hereafter). 

We estimate baryonic cooling rates using similar meth- 
ods to those described by Yoshida et al. (2003), who per- 
formed a detailed hydrodynamic simulation of a 1 h~ 1 Mpc 
volume in order to determine under what conditions gas in 
low mass haloes was able to cool sufficiently to allow col- 
lapse. Our use of a base region almost 500 h _1 Mpc on a 
side allows us to include the important effects of large-scale 
modes (White & Springel 2000; Barkana & Loeb 2004), and 
enables us to model a much earlier generation of baryonic 
cooling in rare haloes with low comoving number density. In 
addition, we model H2 cooling by explicitly integrating the 
time-dependent H2 production rate as described in § 3. 

Population III stars are generally expected to be mas- 
sive (~ lOOM©) due to the large Jeans mass (~ IOOOMq) 
during initial baryonic collapse, and to the inability of zero 
metallicity gas clouds to subfragment into small masses dur- 
ing collapse (e.g. Bromm, Coppi & Larson 1999, 2002; Naka- 
mura & Umemura 2001; Abel et al. 2002; Omukai & Palla 
2003; Bromm & Loeb 2004). We will assume that baryonic 
cooling leads to a single massive star per halo in order to 
make some predictions about their effects on their surround- 
ings and their potential observability. By determining the 
mean separation and clustering properties of the first star- 
forming haloes, we can infer their contribution to reionisa- 
tion. 

If the haloes hosting Pop. Ill stars are highly clustered, 
their associated HII regions could overlap, creating perco- 
lating bubbles of reionisation and, later, metal enrichment. 
Alternatively, if they are weakly clustered, then reionisation 
from the first stars would be confined to isolated HII regions. 
HII regions from Pop. Ill stars have been studied in the red- 
shift range z = 10 — 30 using 1-D codes by e.g. Kitayama et 
al. (2004) and Whalen, Abel, & Norman (2004) who find 
that most UV photons produced by a massive (~ 100M©) 
star are able to escape the virial radius of small haloes such 
as we have simulated here. However, the extent to which 
high redshift ionised regions may overlap with those asso- 
ciated with neighbouring haloes is not yet established. If 
Lyman- Werner radiation or metal enrichment from the first 
generation of stars is effective at suppressing the formation 
of further Pop. Ill stars (e.g. Omukai & Nishi 1999; Glover 
& Brand 2001), overlap will not occur and the original HII 
region will simply collapse once the stars that power it reach 
the end of their lifetime. 

As a first step in distinguishing between these scenar- 
ios, our simulations focus on the largest progenitors at high 
redshifts of a present day cluster region. The structure of 
this paper is as follows: The numerical techniques and as- 
sumptions of these simulations are discussed in § 2. In § 3, we 
describe the criteria required for baryonic collapse. In § 4, we 
show that the first haloes to collapse do so at extremely high 
redshift, and show that after a small redshift interval, many 
other nearby haloes also satisfy the collapse criterion. In § 5, 
we consider the implications of these results; we discuss po- 
tential sources of feedback, and we consider whether this 
early generation of star formation is open to observational 
study. Our conclusions are summarised in § 6. Throughout, 
we assume a flat ACDM model with the following cosmo- 



Table 1. Simulation parameters. Each simulation goes from red- 
shift z s t ar t to redshift z nn . M par t is the particle mass in the high- 
resolution region, M na io the mass of the main halo and r so f t the 
gravitational force softening. 
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logical parameters, which are consistent with the combined 
WMAP/2dFGRS results (Spergel et al. 2003): matter den- 
sity, f2 m = 0.3; dark energy density, = 0.7; baryon den- 
sity, Slbaryon = 0.04; fluctuation amplitude, erg = 0.9; Hubble 
constant h — 0.7 (in units of 100 km s~ Mpc~ ); and no 
tilt (i.e. a primordial spectral index of 1). 



2 NUMERICAL TECHNIQUES 

2.1 The simulations 

We use the parallel gravity solvers Gadget (Springel, 
Yoshida, & White 2001) and Gadget-2 (Springel 2005) to 
resimulate the region surrounding a rich cluster identified in 
a cosmological simulation of a 479 h~ 1 Mpc cube (the VLS 
simulation of Jenkins et al. (2001) and Yoshida et al (2001)). 
This cluster and its environment were part of the sample res- 
imulated and analysed by Navarro et al. (2004) and Gao et 
al. (2004a,b,c). We identified the largest dark matter halo 
at z = 5 in this resimulation, and then resimulated it again 
with resolution improved by an additional factor of 200. We 
repeated this procedure three more times stepping progres- 
sively back in redshift at increasing resolution. 

Table 1 lists all our simulations, giving their initial and 
final redshifts, the mass of an individual particle in the high- 
est resolution region, the mass of the main halo at the final 
time (the original cluster in the case of the VLS simula- 
tion), and the gravitational softening employed. The object 
that we picked at z = 5 did not end up as part of the orig- 
inal cluster selected in the VLS simulation but rather as a 
nearby cluster of smaller mass. This accounts for the differ- 
ence in masses quoted for the final "main halo" in the Rl 
and VLS simulations. The highest resolution resimulation 
ends at redshift 49, has a particle mass of 0.55 h _1 M and 
includes a high resolution subvolume of side 90h _1 kpc. We 
continue to follow the same halo in a resimulation ending 
at redshift 29 with particle mass of 29 h _1 M . Initial con- 
ditions were created using the CMBFAST transfer function 
(Seljak & Zaldarriaga 1996), extrapolated to small scales by 
imposing P(k) tx k -3 . 

2.2 Halo properties 

Our virialized haloes are identified using the spherical over- 
density (SO) algorithm (Lacey & Cole 1994), assuming the 
spherical tophat model (Eke, Cole, & Frenk 1996) in which 
the ACDM virial overdensity, A v i r , in units of the mean den- 
sity is 178 at high redshifts when Q m — 1- For a complex 
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mass distribution such as that expected to arise from the 
flat fluctuation spectrum of our high redshift simulations, 
SO has an advantage over the simpler friends-of- friends halo 
identification algorithm (Davis et al. 1985) in that it is less 
likely spuriously to link together neighbouring haloes or to 
misclassify highly ellipsoidal but unvirialized structures as 
haloes. 

The temperature of gas in virial equilibrium within a 
dark matter halo is given by 
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e.g. Eke et al. (1996), where fi is the mean molecular weight 
(in units of proton mass) and M V i r is the mass within the 
virial radius, r V i r . Simulations by Machacek, Bryan, & Abel 
(2001) and Yoshida et al. (2003) confirm that the local gas 
temperature is relatively close to the virial temperature of 
high redshift haloes, even within the central density peaks. 

Throughout the paper, processes occurring within the 
haloes are calculated assuming uniform density and tem- 
perature within r v i r because they depend in complex ways 
on mixing processes which cannot be treated realistically 
without simulation. This includes the assumption that gas 
is shock-heated upon initial infall. We neglect the possibil- 
ity that some infalling H2 may be destroyed during shock- 
heating, but note that H2 production is rapid enough that it 
would be quickly replenished to a level that allows baryonic 
cooling. Moreover, collisional H2 dissociation rates are small 
below 10 4 K, so any shock heating strong enough to cause 
dissociation would be accompanied by ionisation, providing 
free electrons that would increase the eventual H2 fraction. 
Detailed hydrodynamic simulations are required to model 
the heating of infalling gas accurately. We assume that the 
baryonic mass fraction in the halo is equal to the univer- 
sal mean value. Lower baryon fractions are possible due to 
the suppression of linear baryon fluctuations on small scales 
before the decoupling epoch (e.g. Yamamoto, Sugiyama, 
& Suto 1998; Signh & Ma 2002; Yoshida, Sugiyama, & 
Hernquist 2003). However, even for haloes somewhat below 
the linear theory Jeans mass, gravitational infall into non- 
linear dark matter halo potential wells allows the baryons 
to largely catch up to the dark matter (e.g. Tegmark et 
al. 1997). Yoshida et al. (2003) were able to show that all 
these approximations are reasonably accurate for determin- 
ing whether a halo undergoes sufficient H2 cooling for bary- 
onic collapse of primordial haloes to occur. 

One inevitable uncertainty involves the masses and 
number of the objects formed in the collapse. In order to 
follow the evolution of the baryonic component to the point 
at which individual stars are formed, fully three dimensional 
simulations with radiative transfer would be required. Such 
simulations are currently beyond the limit of numerical sim- 
ulation techniques (Abel et al. 2002; Bromm, Coppi & Lar- 
son 2002; see Bromm & Larson 2004 for a review of recent 
progress) . We thus assume a fiducial collapse model in which 
the baryonic cooling of a halo results in a single ~12OM0 
star. 
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Figure 1. Mass growth of the largest 20 haloes identified in sim- 
ulation R5 at z = 49. The mass plotted is that of the largest 
progenitor identified at each time. 



3 STAR FORMATION CRITERIA 

We apply the methodology of Tegmark et al. (1997) and 
Yoshida et al. (2003) to determine if the baryonic content of 
a halo is able to collapse and form stars. The two require- 
ments are: (1) the cooling rate due to molecular hydrogen 
must be fast enough for the halo to cool within a Hubble 
time; (2) the dynamical heating rate due to mass accretion 
must be slower than the baryonic cooling rate. We assume 
that efficient baryon cooling leads to the formation of an 
increasingly dense gas core that eventually becomes Jeans 
unstable (e.g. Abel, Bryan, & Norman 2000). 

Fig-Qshows that halo growth-rates are remarkably fast, 
the haloes typically doubling in mass over a redshift interval 
Az ~ 2. Here, the growth of the 20 largest haloes in simu- 
lation R5 is illustrated by tracking the evolution of their 
largest progenitor. We show later that because of the rapid 
halo growth, the associated dynamical heating causes a de- 
lay in the baryonic collapse for most haloes. Gas cooling is 
assumed to depend purely on the molecular hydrogen abun- 
dance, which we calculate by integrating the density and 
temperature-dependent H2 production rate over the lifetime 
of each halo, as we now discuss. 

The amount of molecular hydrogen in each dark 
matter halo is computed according to the prescription 
of Tegmark et al. (1997). The primary molecular hy- 
drogen production reaction, H + H - — > H2 + e~ , requires 
the presence of free electrons to form H~. The H2 pro- 
duction rate, dfH 2 /dt ex T 0,88 (where fH 2 = nH,/nH and 
n H = n[H] + n[H+] + 2n[H 2 ] + n[ET] ~ n[H]) is calculated 
as 

— - ~ k m n H x c , (2) 
where k m is the production rate, which is limited mainly by 
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the supply of free electrons for the reaction H + e - — > H~ + 
hv. H2 production depends on the density and temperature 
of the gas, which are taken to be the volume- averaged virial 
quantities, as well as on the free electron fraction x e . We use 
x e = xo = 10 -4 , as obtained for a halo that virializes at z 
~ 50 from the tophat collapse approximation described in 
Tegmark et al. 

We calculate fH 2 using a time-integrated production 
rate. Other authors have taken a different approach, assum- 
ing that the production of H2 stalls at a temperature depen- 
dent value when the remaining free electrons have largely 
recombined (e.g. Tegmark et al. 1997; Yoshida et al. 2003) 
but, as we show in Appendix A, our haloes grow too quickly 
to reach this state. Since halo growth is so rapid, we assume 
that infalling gas mixes uniformly with existing baryons in 
the halo, and this has two primary effects. The first is that 
the H2 within a halo is diluted by infalling material. The 
second is that electrons depleted through recombination are 
quickly replenished by infalling gas such that x c ~ xo. The 
production of H2 in a halo is then described by 
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where < M-^ 2 - > is the average over adjacent simulation 
outputs ti and tj_i, S1h is the hydrogen density in units of 
the critical density, and M is the halo virial mass. The pri- 
mordial H2 fraction, fH 20 is quite small ~ 10 -6 due to H^~ 
photodissociation by CMB photons at high redshift (Anni- 
nos & Norman 1996; Galli & Palla 1998; Nishi & Susa 1999), 
and so can be ignored. The H2 abundance at time ti is thus: 
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We consider a halo to be a star-forming candidate when 
the molecular hydrogen fraction becomes large enough that 
its baryons can cool within a Hubble time: 



7~cool,H 2 



1 k B T 



7 - 1 PA(T) 



< THubblc 



(5) 



where ks is Boltzman's constant and 7 is the adiabatic index 
of the gas. The cooling time, Tcooi,H 2 > is calculated using the 
cooling function, A(T), from Galli & Palla (1998) for molec- 
ular hydrogen rotational line transitions. Fig. [5] compares 
the evolution of fH 2 with the critical abundance, fH 2 ,crit, 
required for cooling in the Hubble time. Sufficient H2 for 
baryonic cooling is produced at z = 50 if we assume di- 
lution by newly-added gas, at which point the halo has a 
mass of 10 5 Mq. However, we will now demonstrate that the 
dynamical heating due to increasing halo mass and the cor- 
responding increasing virial temperature delays collapse. 

Our second criterion for collapse is that the rate of dy- 
namical heating must be lower than the rate of H2 cooling. 
If the dynamical heating rate is larger than the molecular 
hydrogen cooling rate, then nett cooling will not occur and 
baryonic collapse will be prevented. Following Yoshida et 
al. (2003), we require that: 
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(6) 



where Q is the thermal energy per hydrogen atom. The dy- 
namical heating rate is given by 




Figure 2. Evolution of fjj 2 plotted as a function of T v j r (redshift 
decreases from left to right) for the largest halo, and compared 
with fjj 2 criti which is the H2 fraction required for the baryonic 
cooling time to be shorter than the age of the universe. f°H 2 critj 
is shown by solid curves for redshift 70, 60, 50, 40, and 30. To be 
able to cool, the halo must lie above these curves. The solid line 
with data points shows fjj 2 computed by integrating dfjj 2 /dt, and 
assuming that the H2 abundance is diluted by infalling material 
that is both pristine and mixed uniformly into the halo (see Sect. 
3). In all cases, it is assumed that infalling gas is shock- heated 
such that T Eas = T v i r - 
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and the H2 cooling rate is 

dQH 2 ,cool 



1 dt ' 



dt 



A(T)f H2 n H . 



(T) 



(8) 



The quantity dQd yn .hcat/dt is calculated by tracking the 
virial temperature (estimated using Eqn. of the most 
massive progenitor halo through the merger tree. The most 
massive progenitor is defined as the largest halo at the pre- 
vious timestep for which at least half the mass and the most 
bound particle are part of the target halo. In Fig. [3] the 
dynamical heating rate falls below the H2 cooling rate after 
fH 2 ,crit is reached, at a redshift of 47 when the halo mass is 
2.4 x lO 5 h _1 M0. This suggests that baryonic collapse in the 
first halo is delayed by Az = 3, or 4 Myr. Note that there is 
some uncertainty in the redshift of collapse due to the large 
variations of DM/Dz, and also due to the fact that the halo 
mass is not well denned by the SO identification when halo 
shapes are non-spherical. 



4 BARYONIC COOLING IN THE 
SIMULATION 

4.1 The most massive halo 

Summarising our main conclusions from the previous fig- 
ures, the most massive halo in our simulation meets all the 
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Figure 3. Mass growth rate of the largest halo (DM/Dz, solid 
line) compared to the critical growth rate (DM/Dz cr ; t , dashed 
line). DM/Dz cr ; t is defined as the growth rate for which the dy- 
namical heating rate exceeds the H2 cooling rate (Eqn. [fjj. The 
halo growth rate must be below DM/Dz cr i t for the baryons to 
experience nett cooling. The large square (triangle) denotes mea- 
sured (critical) mass growth rate at z = 45. H2 cooling is out- 
paced by dynamical heating in this halo until its mass reaches 
2.4 X 10 5 h- 1 M Q at z = 47. 



criteria for baryon collapse at z ~ 47. At this point, the 
halo mass is 2.4 x lO 5 h _1 M0 and has virial temperature of 
2300K. In this case, dynamical heating of the halo has de- 
layed collapse from redshift 50 to redshift 47, when the uni- 
verse is 51 Myr old. The halo is then expected to undergo 
a runaway collapse that results in the formation of the first 
star or stars. Simulations by Bromm & Loeb (2004) suggest 
that a collapsing primordial star may be able to grow to 
12OM0 in as little as 10 5 yr. Assuming a mass of ~ 120 Mq, 
the first supernova could happen 2.5 x 10 6 years (Schaerer 
2002) after the onset of collapse, by redshift 45. Note that 
the fiducial stellar mass is just below the predicted range 
for pair instability SNe (14O-26OM0), but is still expected 
to undergo an energetic ( Si 10 51 ergs) envelope ejection via 
pulsational pair instability, ultimately forming a black hole 
(Heger & Woosley 2002). 

4.2 Smaller haloes 

In this section, we contrast the cooling of gas in the largest 
halo with subsequent cooling in other large haloes in the sur- 
rounding region. Fig. [I] compares the properties of the 100 
most massive haloes in the simulation at z = 45 with the two 
criteria required for gas to cool and collapse. Although 5 of 
the largest haloes have sufficient H2 for rapid cooling at red- 
shift 45, their dynamical heating rates are still so high that 
only 2 haloes are able to cool by this time. The mass of the 
smaller of these 2 haloes is ~ 10 5 1i -1 Mq. All of these and 8 
other haloes meet both baryonic cooling criteria by redshift 



40, and are thus able finally to cool, as reflected in Fig. 
At redshift 40, the characteristic collapse mass for cooling 
haloes, which we define as the mass above which half of the 
haloes meet both baryonic cooling criteria, is approximately 
2xl0 5 h~ Mq. This characteristic mass increases slowly to ~ 
3xlO 5 h _1 M by z = 29, at which time 80 haloes in the sim- 
ulation are capable of cooling (Fig.[(jJ. Most haloes in Figs. 
14161 lie along a rough path that follows a general sequence 
of increasing mass from lower-left to upper-right covering 
1.5xlO 4 -5xlO 5 h _1 M in Fig.[|and 6x 10 4 -6x 10 7 h _1 M Q 
in Fig. [U Collapse is delayed by dynamical heating for the 
majority of our haloes, as evidenced by the fact that the 
path of the mass sequence crosses the TTj u bble/TH 2 coo i ~ 1 
line before reaching the Qh 2 cool /Qdyn.heat = 1 line. 

4.3 Atomic hydrogen cooling 

Cooling by atomic hydrogen is possible as early as redshift 
36 when the most massive halo reaches a virial temperature 
of 10 4 K, 1.2 x 10 7 years after it formed its first star. At 
this point, if gas has not all been expelled by prior energy 
injection, for example from a supernova, the halo should 
experience much more rapid baryonic cooling. If the high 
efficiency of metal-free atomic cooling leads to efficient star 
formation, this could lead to the formation of an extremely 
early galaxy-like system. The star formation efficiency and 
stellar IMF will depend on the details of the molecular cool- 
ing (or dust and metals if the gas is previously enriched) 
required to cool the gas below 8000K, although isothermal 
contraction could also lead to star formation (Omukai 2001). 
By redshift 29, 13 haloes are hot enough for atomic cooling 
to occur, and a small "group" of galaxies may already have 
formed. Note that this would occur well before the first stars 
form in the simulations by Abel et al. (2002) , Bromm, Coppi 
& Larson (1999, 2002), and Yoshida et al. (2003), but it is 
critically dependent on what happens to the gas after the 
formation and evolution of the first stars. 



4.4 Abundance and spatial distribution of the 
first-star forming haloes 

We cannot use the simulation to determine the abun- 
dance of the first star-forming haloes directly be- 
cause the high-resolution subvolume encompassed only 
90h _1 kpc (comoving) in diameter. The Press & Schechter 
(1974; P-S) prediction for the number density of ob- 
jects more massive than 2.4 x 10 5 /i _1 Mq at z = 47 is 
n P _s(> m)~6x 10" 4 h 3 Mpc" 3 , while the Sheth & Tor- 
men (1999; S-T) formulae give a much larger value of 
5.4 x 10 -2 h 3 Mpc -3 . We do not know a priori which of 
these should be preferred since neither formula has been 
tested against N-body simulations in the regime that is rel- 
evant here. Reed et al. (2003) found the S-T mass function 
to be reasonably accurate even for rare massive objects at z 
= 10, and Springel et al. (2005) showed that the S-T formula 
is much more accurate than the P-S formula out to similar 
redshift. However, these simulations are still orders of mag- 
nitude away from probing the mass range relevant to this 
study. Nevertheless, in what follows we will adopt the S-T 
predictions, but the reader has been warned that this extrap- 
olation is uncertain. Somewhat surprisingly, G05 found that 
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Figure 4. The two criteria for baryonic cooling, 

1"Hubblc/TH 2 cool > 1 (Eqn. |3 and |Qh 2 cool /Qdyn.heatl > 1 

(Eqn.EJ, for the 100 most massive haloes at redshift 45. Haloes 
in the upper right quadrant satisfy both cooling criteria and are 
thus expected to undergo baryonic collapse. Haloes in the upper 
half have enough H2 such that their baryonic cooling time is 
less than the age of the universe. Haloes in the right half have 
H2 cooling rates that are larger than their dynamical heating 
rates and thus experience nett cooling. Dynamical heating delays 
baryonic collapse (haloes in upper left panel) such that only 2 
haloes are capable of cooling by redshift 45, even though 5 haloes 
have sufficiently short H2 cooling times. 
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Figure 5. As Fig. [I] but for z = 40. Now, 13 haloes meet both 
cooling criteria. 



Figure 6. As Figs.l4landl5lbut for the 500 largest haloes present 
at z = 29. Here, 80 haloes meet both cooling criteria. 



the extended Press-Schechter theory (Bower 1991; Bond et 
al. 1991) in the form worked out by Mo & White (1996) 
does accurately reproduce the abundance of haloes in the 
rare high density regions tracked in our R4 and R5 simula- 
tions over the range 30 < z < 50. 

It is interesting to consider how rare is the initial over- 
density that led to the halo where baryons were first able to 
cool in our simulation. To estimate this, we identify the La- 
grangian volume V of the halo by tracing the particles that 
make up the halo at z = 45 (which would correspond to the 
lifespan of a ~ 100M© star formed at z = 47) back to the 
starting redshift, z star t, of the simulation. The Lagrangian 
volume V of such a rare halo is nearly spherical because it 
corresponds to a high peak in the initial density field (e.g. 
Bardeen et al. 1986) . Therefore we can characterise V by its 
mean tophat linear overdensity, 5, and radius, L. We com- 
pute the expected rms density fluctuations on scale L, (j(L), 
from our initial power-spectrum. At z star t = 599, when the 
fluctuation amplitude is very nearly Gaussian-distributed, 
S/a(L) = 6.5, i.e. the halo forms from a 6.5er fluctuation. 

The clustering properties of the first metal-free stars 
will determine their role in the reionisation and metal en- 
richment of the surrounding gas. If these stars are highly 
clustered, then their HII regions could overlap and perco- 
late leading to growing patches of reionisation. If, on the 
other hand, they are more uniformly distributed, percola- 
tion will be delayed, and the eventual reionisation will be 
more uniform. G05 used the Mo & White (1996) formalism 
to estimate spatial correlation functions at z = 47 for haloes 
at least as massive as our largest halo. They found a bias 
factor of 24 which results in a comoving correlation length of 
2.5 h _1 Mpc, only slightly smaller than that of star-forming 
galaxies today. The comoving abundances estimated above 
are lower than those of such galaxies, however, and the lu- 
minosities predicted for our first star systems are very much 
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lower. As a result, the first light sources are relatively weakly 
clustered despite their high bias. Neighbouring haloes that 
lie close enough to allow the possibility of overlapping HII 
regions are extremely rare at redshift z = 47. As we discuss 
next, the HII regions of haloes in the region we have followed 
may, nevertheless, percolate. 



5 DISCUSSION 

We have shown that gas is able to cool and collapse in haloes 
of mass 2.4x lO J h _1 M0 already at z = 47, and that this is 
followed quickly by the collapse of the gas in less massive 
haloes in the surrounding region. In just another 22 Myr, 
at z = 36, the first haloes able to cool by atomic cooling 
have formed. In this section, we consider the impact of the 
collapse of the first haloes on their surroundings, focusing in 
particular on the size of the ionisation patches that are likely 
to be formed. We consider the observable consequences of 
the formation of the first stars in these haloes at the end of 
this section. 

5.1 Ionising radiation from the first stars 

We begin by calculating the role of the first stars in ionising 
the surrounding intergalactic medium (IGM). In the next 
subsection, we will consider the effect of the first stars on 
the star formation rate in surrounding haloes. In order to 
estimate limits on the impact of this first baryonic collapse, 
we use the flux and lifespan calculations for Pop. Ill stars 
by Schaerer (2002). We first establish some fiducial values. 
The efficiency of production of ionising photons peaks at a 
stellar mass of 12OM0. We parameterise the total ionising 
photon flux resulting from the collapse of the gas in a halo 
as 

n io niso = 1.1 x 10 64 mi2ofi20, (9) 

where mi2o is the mass, in units of 12OM0, of ionising stars 
formed per halo, and fi2o is the total number of ionising 
photons emitted per stellar mass relative to the specific ion- 
ising flux from a 12OM0 star. Since fi2o is a weak func- 
tion of stellar mass, ni on j sc is a weak function of the ac- 
tual mass of individual star(s) formed, and is mainly de- 
pendent on the total mass turned into stars. Later, we will 
assume that the 2.4 x 10 h _1 M0 dark matter halo (which 
contains 3.2 x 1O 4 M0 baryons), forms a single I2OM0 star 
(i.e. mi2ofi20 = 1). The maximum possible size of the HII 
region around such a star-forming object is 76(mi2ofi2o) 1//a 
kpc (comoving), if we assume no recombination (i.e. the fi- 
nal number of ionised atoms is equal to the total number of 
ionising photons produced by the star) and that the sphere 
expands into the mean baryon density. In practise, the high 
densities at such high redshift give rise to high recombina- 
tion rates that are likely to limit the HII region significantly. 

The size of the ionised region is uncertain for several 
reasons. Firstly, the recombination rate depends sensitively 
on the local gas density. For simplicity, we assume that gas 
traces dark matter, which is likely to be a good approxi- 
mation outside the virial radius in the absence of feedback 
effects, although there is little reason to expect this to be 
correct within the virial radius. Secondly, a further com- 
plication arises from the effect of Compton cooling on the 



gas. The Compton cooling timescale is shorter than the re- 
combination timescale in low density regions at high red- 
shifts. This could serve to quench the growth of ionised re- 
gions because the recombination rate increases as gas cools. 
Thirdly, the expected short lifetimes of these massive stars 
could result in only a fraction of haloes hosting their first 
star at any one time. However, after the first star's death, 
the enhanced electron fraction within its HII region may 
promote subsequent star formation, albeit of stars possibly 
of lower mass (O'Shea et al. 2005). Fourthly, supernovae 
from the first stars may create bubbles larger than their 
HII regions. Finally, as we discuss in § 5.2, Lyman- Werner 
radiation may destroy molecular hydrogen in neighbouring 
haloes, thus limiting the number of stars that form in H2- 
cooled hosts, although stars in atomic-cooled haloes may 
soon dominate the stellar population. In practice, of course, 
the ionised region will not be spherical but will depend on 
the morphology of the surrounding gas distribution. 

We have estimated the size of an HII region using 
a steady-state 1-D code. The recombination rate is gov- 
erned by the rate coefficient, k rccomb = 1.88 x io 10 T" a64 
(cm 3 s _1 ) (Hutchins 1976), and the temperature of the 
ionised gas is set to T = 10 4 K. The size of the HII region is 
taken simply to be the radius at which the recombination 
rate within the sphere is equal to the ionising photon pro- 
duction rate. We assume that the gas follows an NFW profile 
(Navarro, Frenk & White 1996, 1997) with a concentration 
of 1.75 at z = 47, and 1.87 at z = 29, which provides a good 
fit to the dark matter profile in the simulation (G05). To 
estimate the effects of the shock front, we smooth the gas 
in the central portion of the HII region to have a uniform 
density over the sound crossing distance for a star with a 2.5 
Myr lifespan (assuming ionised gas at 10 4 K). This smooth- 
ing also makes the estimated HII region size less sensitive 
to the assumed profile shape. The mass density surrounding 
the largest halo in the simulation approximately follows an 
NFW profile extrapolation until it reaches ~ 10 times the 
mean density at ~ 3r v i r , beyond which we assume a uniform 
density. 

For the fiducial mi2ofi20 = 1 case, at redshift 47, the HII 
region extends just beyond the virial radius, to 1.2r v ; r in our 
calculation. However, by redshift 29, physical densities are 
low enough that over half the ionising photons from a single 
I2OM0 star will escape a halo of the characteristic mass for 
H 2 -cooled haloes at this redshift (~ 3 x lO 5 h _1 M ). Addi- 
tionally, a much larger number of haloes are now susceptible 
to H2 cooling, raising the possibility that percolation of HII 
regions could create a large HII bubble. At z = 29, the HII 
region around a single 12OM0 star would extend to 9r v i r 
(21 h _1 kpc, comoving). This is large enough for the entire 
~ 75 h _1 kpc (comoving) radius region that we have simu- 
lated to have an HII filling factor of ~ 1, provided that each 
halo maintains a similar sized HII region. Of course, such 
regions are rare and the contribution to the ionised volume 
of the universe as a whole is small (see § 5.3). The percola- 
tion of these HII bubbles is illustrated in Fig. |7| left panel. 
This scenario is probably unrealistic due to short recombi- 
nation times of ~ 1 Myr in the high density ( £ 10 times 
the mean density) vicinity of these haloes, and due to po- 
tential removal of gas by feedback. The short recombination 
times imply that percolation of ionised bubbles will be sig- 
nificantly reduced if a new ionising source(s) does not form 
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before the end of the first star's lifetime. The right-hand 
panel in Fig. [7| illustrates how percolation will be reduced if 
each HII region survives only one recombination time after 
the lifetime of the primary star. Finally, we have considered 
whether Compton cooling could affect the size of individual 
HII regions by enhancing the recombination rate, and find 
that this effect is likely to be unimportant over the lifetime 
of a massive star. At redshift 29, the 1.4 Myr Compton cool- 
ing timescale is of similar order, but generally longer than 
the typical local recombination time. 

Thus, while the size of an HII region is quite uncer- 
tain, it seems possible that a large bubble around the first 
star-forming halo could become ionised well before the rest 
of the universe. One might then expect similarly large and 
widely separated HII bubbles to be scattered throughout the 
Universe prior to reionisation. 



5.2 Feedback from the first generation of stars 

In § 4, we neglected the impact of the collapse of the first 
stars on neighbouring haloes. The picture we have presented 
so far only takes into account the ability of a halo to cool 
and form stars according to local criteria without taking into 
account external influences. We now briefly discuss how neg- 
ative feedback from the first stars may impact star formation 
in the surrounding region. 

There are three (possibly more) feedback mechanisms 
that could mean that the ability of a halo to cool is not a 
local process, and does not depend solely on the mass of the 
dark matter halo and its formation history. 

• Lyman-Werner radiation: Radiation in the Lyman- 
Werner bands can dissociate hydrogen molecules. If the 
destruction rate exceeds the rate of production, the halo 
will not build up the amount of H2 needed to cool. As 
a result, the formation of the first star may "sterilise" 
the surrounding region and prevent other stars from form- 
ing (e.g. Haiman, Rees, & Loeb 1997; Ciardi, Ferrara, & 
Abel 2000; Haiman, Abel, & Rees 2000; Ricotti, Gnedin, & 
Shull 2002; Mackey, Bromm, & Hernquist 2003; Yoshida et 
al. 2003). 

The photo-dissociation timescale for H2 is given by 



t — V- 1 F" 1 

^diss K diss r trans 

k diss = 1.1 X 10 8 j L w(s _1 ) 



(10) 

(11) 



(Abel et al. 1997), where ji,w is the Lyman-Werner flux at 
hv = 12.87eV in units of 10" 21 ergs s _1 cm" 2 Hz" 1 . In the 
absence of self-shielding, the transmission factor F tr ans is 
equal to 1. 

It is convenient to parameterise the dissociation timescale 
in terms of the photo-dissociating luminosity (12.24ev- 
13.51ev) for the fiducial 120M o star L^: 



tdis 



1.3 



1 f d 



L^w ^lMpc 



(it;)' 



Ft, 



(Myr), (12) 



where d is the distance from the star in comoving Mpc. 

In the absence of self- shielding (see below) , a single 120Mq 
star produces enough Lyman-Werner flux to dissociate H2 
within a volume of ~ 1 comoving Mpc 3 . (Here, the "steril- 
isation by dissociation" radius is determined by estimating 



the distance at which the dissociation timescale is compara- 
ble to the H2 formation timescale, which is typically ^ 1 
Myr.) 

This calculation ignores the fact that the surrounding 
haloes are likely to contain dense cores. If the optical depth 
in these regions is sufficiently high, the inner mass will es- 
cape dissociation because of self-shielding in the outer parts 
(e.g. Glover & Brand 2001; Kitayama et al. 2001; Machacek 
et al. 2001). Unfortunately, the effectiveness of self-shielding 
depends on the internal velocity structure of the halo gas 
caused by infall and turbulent motions. The properties of 
these flows are uncertain and they will be difficult to model 
adequately even in the highest resolution simulations. 

If the gas velocities are negligible, the transmission factor 
is given approximately by: 



Ft,. 



1, 



Nt 



10 14 cm" 2 



-3/4' 



(13) 



(Draine & Bertoldi 1996), where Nh 2 is the H2 column 
density, integrated from outside-in. Because it neglects gas 
motions, this calculation gives the maximum possible self- 
shielding (Machecek et al. 2001). For an NFW halo pro- 
file with a concentration of 1.75, and assuming that the H2 
density traces the mass density, Eqn. 1131 implies that 90% 
of the H2 mass (r > 0.2r v i r ) will be prone to dissociation 
by a single 12OM0 star located at a distance of 65 kpc at 
redshift 45, for a typical H2 fraction of 10 -4 . This implies 
that dissociation is still important at large distances. For the 
typical halo separation in our simulations (~ 10 h~ 1 kpc), a 
transmission factor of less than 1.5xl0~ 4 would be required 
to ensure nett H2 production (tdiss ~ lMyr) at redshift 45. 
However, at this distance, the central baryonic mass that 
is adequately shielded is much less than that required to 
form a massive star unless either the baryon density or the 
H2 fraction are significantly enhanced. On the other hand, 
baryonic cooling is likely to concentrate baryons strongly to- 
wards the halo centre (e.g. Yoshida et al. 2003), increasing 
the Lyman-Werner optical depth and decreasing the effect of 
Lyman-Werner feedback. These issues can only be addressed 
with detailed numerical hydrodynamic simulations. 

We estimate that neither L-W shielding by intergalactic 
H2 nor by intergalactic neutral Hydrogen will be important. 
Using the transmission factor of Eqn. 1131 and assuming a 
primordial H2 fraction of 10 -6 at redshift 50, L-W shielding 
by intergalactic H2 will be negligible except over very large 
distances (F tra ns = 0.1 at 1 Mpc, comoving). This estimate 
likely over-predicts the inter-galactic shielding effect because 
doppler shifts due to Hubble expansion and intrinsic gas mo- 
tions will decrease the optical depth. More importantly, any 
shielding by inter-galactic H2 will be short-lived because it 
will quickly be dissociated. Shielding by inter-galactic neu- 
tral hydrogen is also ineffective because L-W absorption will 
be largely confined to the closely spaced lyman series lines 
in the upper energy L-W bands. Even at a distance of 1 
Mpc (comoving), less than half of total L-W radiation will 
be absorbed by neutral H at redshift 50. 

We have rerun our H2 production model including the 
dissociating effects of Lyman-Werner radiation under the as- 
sumption that self-shielding is unimportant. The L-W flux 
from a single 120 Mq star is enough to quickly sterilise the 
entire (75 h _1 kpc comoving radius) high resolution subvol- 
ume of the R4 simulation. However, if we then switch off the 
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a) HH region size is maintained b) HII can recombine after 
by subsequent source(s) 1 st star switches off 




Figure 7. In this schematic image, an HII region powered by a 120 Mq Pop. Ill star is centred on each z = 29 R4 halo that meets the 
H2 cooling criteria of § 3 at redshifts 45, 40, and 29. In the left panel, the 21 kpc (comoving) radius HII regions (see text) are assumed 
to be maintained after the end of the first star in each halo by subsequent ionising sources. In the right panel, HII regions are assumed 
to survive for only one recombination time (~ 1 Myr) beyond the life of the first star in each host. Significant local percolation of HII 
regions is seen in both cases. In § 5.2, we discuss feedback mechanisms that may greatly alter the number of star-forming haloes in 
regions such as this. 
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L-W flux after a 2.5 Myr stellar lifetime, a second halo is able 
to produce enough H2 to meet the baryonic cooling criteria 
within 1-2 Myr. Thus, the sterilising effects of a single mas- 
sive star are temporary, but may nevertheless limit the star 
formation density to ~ 1 star per comoving Mpc per Myr 
in H2-cooled hosts. However, if self-shielding within dense 
baryonic cores can cause a significant time delay between 
the onset of baryonic collapse and subsequent sterilisation 
of neighbouring haloes, then there may exist a brief win- 
dow during which multiple haloes will have time to produce 
enough H2 for baryonic cooling. A sufficiently long window 
could allow a rapid "mini-burst" of H2 induced population 
III star formation in a previously sterilised region. 

Thus, depending on the detailed structure of collapsing 
clouds, the radius out to which Lyman- Werner feedback is 
effective could be as large as ~ 1 Mpc, which would sub- 
stantially limit the number of haloes that undergo collapse 
via H2 cooling. Self-shielding is likely to mediate this feed- 
back to some degree, but it is probably unimportant except 
for dense baryonic cores. However, we stress that even if all 
the H2 in a halo were destroyed by external Lyman-Werner 
radiation, atomic cooling alone should still be sufficient to 
trigger baryonic collapse in our haloes once they reach a 
sufficiently high virial temperature (e.g. Omukai 2001). 

• X-ray and far UV radiation: The first supernovae and 
mini-qso's will bathe the surrounding region with ionising 
photons (e.g. Madau et al. 2004). Those of sufficiently high 
energy have small ionisation cross-sections and will escape 
from the HII region. However, even if they only ionise a 
small fraction of the neutral hydrogen in a nearby halo, the 
increase in the free electron density will catalyse the produc- 
tion of H2. If this process is significant in the centres of these 
haloes, it would promote the collapse of their gas, acting as 
a positive feedback mechanism (e.g. Haiman, Rees, & Loeb 
1996; Haiman, Abel, & Rees 2000; see however, Machacek, 
Bryan, & Abel 2003, who find that a uniform X-ray back- 
ground has little effect on H2-induced cooling and collapse 
of haloes in simulations) . 

• Metal enrichment and SNe: The energy released by the 
explosion of the first SNe sends a blast wave into the sur- 
rounding IGM. If cooling haloes are fragile, the associated 
mechanical energy can remove portions of their outer enve- 
lope and prevent further cooling (e.g. Rees 1985; Bromm, 
Yoshida, & Hernquist 2003; Kitayama & Yoshida 2005). 

Pop. Ill stars with masses in a narrow range around 
2OOM0 are thought to end their life as a pair production 
supernova (e.g. Heger & Woosley 2002). In this case, the su- 
pernova ejecta is very metal rich and, if this material mixes 
in with the gas in nearby haloes, the overall gas cooling ef- 
ficiency will be enhanced. Metal enrichment is, in any case, 
required in order to enhance baryon cooling at low temper- 
atures (<10 4 K) so as to promote the fragmentation that 
results in the low-mass IMF typical of population II stars 
(Schneider et al. 2002, 2003; Mackey, Bromm, & Hernquist 
2003). At redshift ~ 30, energetic (~ 10 53 ) SNe blastwaves 
may expand up to ~ 20 kpc (comoving) before the initial 
fast expansion is "stalled" (estimated according to Yoshida, 
Bromm, & Hernquist 2004), suggesting that metal enrich- 
ment from neighbouring haloes is possible. 

Thus, supernovae blast waves can produce both positive 
and negative feedback. However, the region affected by these 



processes is likely to be small compared to the region poten- 
tially affected by the Lyman-Werner radiation. 

• H~ Photo-detachment Another possible mechanism 
suppressing H2 formation is the photo-detachment of H~ 
ions before they can combine with a hydrogen atom to 
produce H2. We estimate the radius over which photo- 
detachment suppresses H2 formation by equating the photo- 
detachment rate to the rate at which H~ reacts to form H2 
(H + H~ — > H2 + e~; rate taken from Hirasawa 1969), as- 
suming virial overdensity. We assume the effective tempera- 
ture and luminosity of a 120 Mq Pop. Ill star from Schaerer 
(2002), and employ the Tegmark et al. (1997) fitting function 
to photo-detachment cross-sections from Wishart (1979). 
This yields a photo-detachment distance of 2.5 kpc (comov- 
ing) at z = 29, roughly equal to the virial radius of a typi- 
cal halo on the verge of baryonic cooling at this time, and 
much smaller than the HII region expected from such a mas- 
sive star. Thus we do not expect H - photo-detachment by 
~ IOOMq Pop. Ill stars to impact baryon cooling in neigh- 
bouring halos. 

In summary, our calculations show that Lyman-Werner 
radiation could "sterilise" a region as large as ~ 1 Mpc (co- 
moving) around a primordial massive (120 Mq) star. How- 
ever, the effects of self-shielding are so uncertain that we can- 
not rule out the possibility that dissociation feedback may be 
ineffective at suppressing further star formation. The uncer- 
tainties in other feedback mechanisms are also considerable. 
In the preceding section we discussed the possibility that the 
HII regions around the first stars might percolate giving rise 
to large multi-halo HII bubbles. Whether or not this hap- 
pens depends primarily on how effective the Lyman-Werner 
sterilisation is at suppressing star formation in the vicinity of 
a star, and on the impact of the other feedback schemes dis- 
cussed above. However, even if sterilisation does occur, this 
can only delay star formation in the region, since atomic 
cooling, already possible in 13 haloes in our simulation by 
z = 29, will ultimately become the dominant path for star 
formation. For the case of maximally effective sterilisation, 
the S-T number density of haloes capable of atomic cooling 
surpasses that of H2-cooled haloes (~ 1 per h _3 Mpc 3 ) by z 
~ 25. Gas in atomic-cooled haloes could dominate the to- 
tal stellar mass and UV radiation considerably earlier than 
this if, as seems likely, this gas turns into stars more effi- 
ciently than gas within H2-cooled haloes due to faster rate 
and greater temperature sensitivity of atomic cooling. 

Unfortunately, it is not yet possible to rule out any par- 
ticular scenario. We proceed by exploring the extreme case 
where dissociative feedback is unimportant, and percolation 
leads to the formation of a large ionised region. 

5.3 Contribution to reionisation 

The sequence of events following the formation of the first 
generation of Pop. Ill stars (z £ 45) is uncertain because 
of our poor understanding of the interplay between negative 
feedback from Lyman-Werner radiation and positive feed- 
back from the X-ray/UV background. If these processes are 
unimportant or if positive feedback dominates, our calcula- 
tions suggest that large (~ 100 h _1 kpc comoving) patches 
of the universe will be reionised at z > 30. The first regions 
of the universe to collapse may thus become reionised much 
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earlier than the z ~ 20 inferred from the WMAP microwave 
background radiation data (Kogut et al. 2003). This could 
have dramatic effects on low mass galaxies forming in these 
regions whose gas could be removed suppressing subsequent 
star formation. 

Observational constraints derived from observations of 
high redshift QSOs (Fan et al. 2001) and from the temper- 
ature of the high redshift IGM (Theuns et al. 2002) imply 
that complete re-ionisation occurs at z = 6 — 9. In contrast, 
the high optical depth for Thomson scattering inferred from 
the WMAP data suggests a much earlier epoch for reionisa- 
tion, at least over significant fractions of the Universe. It is 
interesting to consider whether large patches of early reion- 
isation could reconcile these contradictory constraints (see 
e.g. Bruscoli, Ferrara, & Scannapieco 2002; Cen 2003; Cia- 
rdi, Ferrara, & White 2003; Haiman & Holder 2003; Sokasian 
et al. 2004). However, according to the number densities in- 
ferred from the S-T formula, regions such as the one we 
have simulated are quite rare and so they can only make a 
small contribution to global reionisation and to the Thomp- 
son optical depth. Moreover, stars in H2-cooled haloes could 
even delay reionisation if their sterilising effects prevent H2- 
cooled star formation over large distances. Alternatively, 
strong Lyman- Werner feedback from the first population III 
stars might actually expedite reionisation by preventing the 
expulsion of gas by feedback associated with star formation 
in early low mass haloes, leaving plentiful gas supplies when 
haloes grow large enough for efficient atomic cooling. 

5.4 Observability 

Finally, we consider whether the first generation (z £ 45) 
of population III stars that is predicted by our simulations 
could actually be observed. We focus first on the case where 
H2 is unaffected by feedback processes, which are highly un- 
certain. The supernovae resulting from the collapse of these 
stars are potentially observable, particularly if they can be 
detected as gamma-ray bursts (GBRs). The time dilation of 
bursts occurring at higher redshifts implies that the source 
can be observed closer to maximum light, thus compensat- 
ing for the increasing luminosity distance. In a recent paper, 
Gou et al. (2004) showed that GRBs can be detected in X- 
rays out to z = 30 and beyond in exposure times as short as 
300 sec. 

The rest-frame ultraviolet emission from the burst will 
be strongly absorbed shortwards of Ly-a, so that detection 
will need to be made in the mid infrared. Gou et al show 
that the James Webb Space Telescope will have the required 
sensitivity to detect GRBs out to z = 37 in the 4.8/^m band 
in only 1 hour. Detection of GRBs at higher redshift should 
also be possible in longer wavelength bands and larger ex- 
posure times. Radio afterglows may also possible at z ~ 30 
(Ioka & Meszaros 2005). 

The number of haloes massive enough to support star 
formation lying in our past lightcone at z > 45 is of the order 
of 10 9 , according to the S-T prediction. If we assume that 
each of these haloes forms a single massive star with a 2.5 
Myr lifespan whose endpoint is a GRB, then the event rate 
over the entire sky would be of the order of one per month. 
Our estimated abundances (which, as we discussed in § 4.4, 
are uncertain) are higher than, but in rough agreement with 
those of Miralda-Escude (2003) who, on the basis of the 



P-S formula, predicts the highest redshift star in our past 
lightcone to be at redshift 48, assuming that such stars form 
in haloes with a virial temperature of 2000K. Our z ~ 50 
supernovae rates are of similar order to those predicted by 
Wise & Abel (2004) on the basis that supernovae are hosted 
by haloes of mass ~ 10 Mq whose abundances is given by 
the S-T model. Note that we have ignored GRB beaming 
in this calculation; the number of observable GRBs that 
are beamed toward us is likely to be lower by a factor ~ 
100-1000 (e.g. Frail et al. 2001). 

The rate of potentially observable GRBs or SNe forming 
in haloes where gas has cooled via H2 increases rapidly at 
lower redshifts. At z = 29 in our simulation we expect ~ 10 
supernovae per day over the entire sky, or 1 per ~ 10 square 
degrees per year (again, ignoring beaming effects). If sterili- 
sation by Lyman- Werner radiation strongly suppresses star 
formation, the number of observable GRBs at z>29 could 
be reduced by a significant factor, depending on the uncer- 
tain effects of self-shielding. However, even in that case, the 
number density of haloes capable of sustaining atomic cool- 
ing would soon surpass that of H2-cooled haloes, leaving 
open the possibility of abundant pre-reionisation observa- 
tional targets. The increase in the expected event rate at 
lower redshifts suggests that GRBs or SNe should indeed be 
observable in the pre-reionised universe. 

We conclude that if the highest redshift Pop. Ill stars 
end their life as GRBs, the prospects for their detection are 
promising. Even if these stars do not have a GRB phase or 
if the beaming angle is small, observational evidence of the 
first stars may still be available, for example, through direct 
detection of SNe (Panagia 2003) or of emission from an early 
generation of accreting black holes (e.g. Madau et al. 2004). 
The difficulty, however, will lie in identifying such rare and 
faint objects and distinguishing them from foreground con- 
taminants. 



6 CONCLUSIONS 

We have modelled the cooling of gas by molecular hydrogen 
line emission in a series of N-body simulations in order to 
identify the first generation of dark matter haloes capable 
of forming stars. We followed the formation of structure in a 
region selected to host a rich cluster today, repeatedly resim- 
ulating, at increasing resolution, the largest progenitors of a 
massive halo at the final time. To determine whether the gas 
is able to cool, we estimated the production rate of H2 and 
considered the balance between radiative cooling and dy- 
namical heating. The dramatic growth of the halo at early 
times has the nett effect of delaying, but not preventing the 
cooling of gas. Our main results are as follows. 

• The first star-forming haloes: Since we have fol- 
lowed evolution in a special region of the universe (selected 
to have a rich cluster today), we find that formation of the 
first stars occurs significantly earlier than commonly found 
in simulations of volumes too small to sample such rare over- 
densities. At a redshift of 47, there is already a population 
of haloes in the simulation that are undergoing baryonic col- 
lapse via H2 cooling and can host the first stars. 

The mass of these haloes is typically 2.5xlO 5 h _1 M0 
and their comoving space density is comparable to that of 
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dwarf galaxies today. By redshift 29, 80 haloes in the high- 
resolution simulation subvolume of radius 75h _1 kpc are ca- 
pable of forming stars. The largest halo reaches a virial tem- 
perature of 10 4 K, at which gas in equilibrium can begin to 
cool by more efficient atomic processes, at z = 36, poten- 
tially giving rise to a first generation galaxy. By z = 29, 13 
haloes in the simulation have grown large enough to sustain 
atomic cooling. 

• Ionising effects: The formation of the first stars will 
ionise the gas around them. At z = 47, the hydrogen recom- 
bination rates are high enough that these Hll regions are 
likely to be contained almost entirely within their own host 
halo. However, by z = 29, gas densities have dropped suffi- 
ciently that massive (~ 100 Mq), population III stars are 
able to ionise regions significantly beyond the halo. Perco- 
lation of neighbouring regions could lead to the formation 
of an HII superbubble. However, such regions would be too 
widely separated to contribute significantly to reionisation 
at this time. 

• Feedback: Feedback effects that could alter the cool- 
ing properties of the gas are likely to occur but their nature 
(and, in some cases, even their sign) is difficult to calcu- 
late. Lyman- Werner radiation from just a single massive 
star could sterilise a large surrounding volume by dissoci- 
ating H2, severely limiting star formation throughout the 
volume. However, the effectiveness of this feedback chan- 
nel would be limited if the H2 within dense halo cores is 
self-shielded by the outer halo layers. Supernova feedback 
could hinder star formation by removing gas, or could en- 
hance star formation by blast wave compression, by metal 
enrichment, or through high energy radiation which might 
produce electrons that catalyse H2 formation. 

Even if feedback is effective enough to suppress H2 cooling 
within a region of size ~1 Mpc (comoving) around the first 
star, further star formation is still expected in haloes where 
gas can cool by atomic processes. This happens as early as 
z = 36 in our simulated region. 

• Observability: The prospects for observing the very 
earliest generation of stars are encouraging. The James 
Webb Space Telescope may be able to view directly pre- 
reionisation supernovae and gamma-ray bursts from popu- 
lation III objects in the rest-frame mid-infrared (Panagia 
2003). X-ray emission or radio afterglow from these early 
gamma-ray bursts may also detectable (Gou et al. 2004). 
The event rate of supernovae at z = 29 and higher could 
be as large as ~ 10 per day. However, if feedback due to 
Lyman- Werner radiation is important, the event rate could 
be substantially lower. 
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APPENDIX A: 

In this appendix, we discuss the production of H2, compar- 
ing our method for calculating fH 2 with that employed by 
Tegmark et al. (1997) and Yoshida et al. (2004). As the elec- 
trons left over from recombination are depleted, the produc- 
tion of H2 slows down. This leads to the "asymptotic" H2 
fraction, fH 2 ,asym oc T 1 ' 52 , where T is assumed to be equal to 
the virial temperature of the halo (see Eqn. 17 of Tegmark 
et al. 1997). 

Since T vir oc M 2/3 , where M is the halo mass, it fol- 
lows that dfH 2 ,asym/dt 6c dM/dt as the halo grows in mass, 
where dfH 2 ,asym/dt is the H2 production rate needed to sus- 
tain fH 2 ,asym, which increases as the halo grows. The asymp- 
totic H2 abundance is maintained only if the halo growth is 
sufficiently slow so that the actual production rate 

df H2 df H2 ,asym , a1 n 

"dT > dt ' (A1) 

where dfH 2 /dt is the halo H2 production rate. If the halo 
grows too rapidly, molecular hydrogen production will not 
be able to keep up (even with the primordial abundance of 
free electrons) and the cooling efficiency of the halo will be 
reduced. In Fig. lAll we show that the H2 production rate for 
our largest halo drops below the rate required to maintain 
the asymptotic H2 abundance at z = 62. This means that for 
z < 62, fH 2 for the halo increases more slowly than fH 2 ,asym 
because of the high halo growth rate. Similar curves are 
readily produced for other high redshift haloes. It is thus 
necessary to estimate H2 abundances by integrating dfH 2 /dt 
over time, rather than by simply adopting "asymptotic" H2 
abundances (we also include dilution effects due to mergers 
as described in § 3). Some haloes have periods where dfH 2 /dt 
(calculated from Eqn.|2J is larger than dfH 2 ,asym/dt, in which 
case we assume that electron depletion limits H2 production 

tO df H2 ,asym/dt. 
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Figure Al. Comparison of the estimated production rate of 
molecular hydrogen dfjj 2 /dt in a halo (solid curve) versus 
[dfjj 2 / dt] asym ( dashed curve) , the rate of production needed to 
maintain the "asymptotic" H2 fraction. fjj 2 aBym is the result- 
ing H2 abundance after H2 production slows dramatically due to 
the depletion of catalysing free electrons. The growth of the halo 
is too rapid to maintain [dfji 2 /dt] a sym below z ~ 62, and so at 
z < 62, the Ho abundance of the halo will be below ftr., • In 
this regime, the H2 abundance is thus limited primarily by the 
halo's rapid growth rather than by electron depletion. We thus 
determine the H2 abundance from a time integration of the H2 
production rate (see text). The curves that stop at redshift 49 
follow simulation 'R5' until its completion, while the curves that 
continue to redshift 29 are for the same halo, which is modelled 
at lower resolution in simulation 'R4'. 
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